Hyades Member K2-136c: The Smallest Planet in an Open Cluster with a Precisely Measured Mass

Photometric observations of the K2-136 system were collected with the Kepler spacecraft (Borucki et al. 2008) through the K2 mission during Campaign 13 (2017 March 8 to May 27). K2 collected long-cadence observations of this system every 29.4 minutes.

Photometric observations of the K2-136 system were also collected with the Transiting Exoplanet Survey Satellite (TESS) spacecraft (Ricker et al. 2015) during Sector 43 (2021 September 16–October 10) and Sector 44 (2021 October 12–November 6). TESS collected long-cadence full-frame image observations of this system every 10 minutes in Sector 43 and short-cadence observations every 20 s in Sector 44.

We collected 93 RV observations using the High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N) spectrograph (Cosentino et al. 2012, 2014) on the Telescopio Nazionale Galileo (TNG). The first 88 spectra were collected between 2018 August 11 and 2019 February 7 (programs A37TAC_24 and A38TAC_27, PI: Mayo), and the final 5 spectra were collected between 2020 September 18 and October 31 by the HARPS-N Guaranteed Time Observation program. RVs and additional stellar activity indices were extracted using a K6 stellar mask and version 2.2.8 of the Data Reduction Software (DRS) adapted from the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) pipeline. The spectra had an average exposure time of 1776.5 s, and the average signal-to-noise ratio (S/N) on the order around 550 nm was 51.1. The RV standard deviation was 6.9 m s⁻¹, and the RV median uncertainty was 1.6 m s⁻¹. The stellar activity indices extracted and reported in this paper include the cross-correlation function (CCF) bisector span inverse slope (hereafter BIS), the CCF full width at half maximum (FWHM), and S_HK (which measures the chromospheric activity via core emission in the Ca ii H and K absorption lines). The observation dates, velocities, and activity indices are provided in Table 2.

We collected 22 RV observations using the ESPRESSO spectrograph (Pepe et al. 2021) on the Very Large Telescope (VLT) between 2019 November 1 and 2020 February 27 (program 0104.C-0837(A), PI: Malavolta). RVs and additional stellar activity indices were extracted using a K6 stellar mask and the same pipeline as the HARPS-N observations (DRS version 2.2.8). The typical exposure time for the spectra was 1800 s, and the average S/N at order 111 (central wavelength = 551 nm) was 79.9. The RV standard deviation was 7.8 m s⁻¹, and the RV median uncertainty was 0.70 m s⁻¹. These observations and indices are also provided in Table 2.

Near-ultraviolet (NUV) observations of K2-136 were taken as part of a broader Hubble Space Telescope (HST) program observing the Hyades (GO-15091, PI: Agüeros). The target was exposed for 1166.88 s on 2019 September 13 using the photon-counting Cosmic Origins Spectrograph (COS; Green et al. 2012) in the G230L filter and had no data quality flags.

After initial data reduction through the CALCOS pipeline version 3.3.10, we additionally confirmed that the star was not flaring during observations by integrating the background-subtracted flux by wavelength over 1 and 10 s time intervals in the time-tagged data. No flares above 3σ were identified.

K2-136 was the target of an XMM-Newton (XMM) 43 ks observation on 2018 September 11 (Obs. ID: 0824850201, PI: Wheatley). The observation was processed using the standard Pipeline Processing System (PPS, version 17.56_20190403_1200;Pipeline sequence ID: 147121). The source detection corresponding to K2-136 was detected by both the pn and metal oxide semiconductor (MOS) cameras, for a total of 800 source counts in the 0.2–12.0 keV energy band. The X-ray source has a data quality flag SUM_FLAG=0 (i.e., good quality). No variability or pileup were detected for this X-ray source.

One parameter of special interest is the stellar rotation period, which we include as a parameter in our RV model (see Section 4.5). M18 conducted a Lomb–Scargle periodogram on the K2 light curve and reported a rotation period of 15.04 ± 1.01 days. Ciardi et al. (2018) analyzed the same light curve and found a rotation period of 15.2 ± 0.2 days through a Lomb–Scargle periodogram and 13.8 ± 1.0 days through an autocorrelation function. Livingston et al. (2018) conducted a Gaussian process (GP) regression, a Lomb–Scargle periodogram, and an autocorrelation function on the light curve and found a corresponding rotation period of ${13.5}_{-0.4}^{+0.7}$ days, ${15.1}_{-1.2}^{+1.3}$ days, and ${13.6}_{-1.5}^{+2.2}$ days, respectively. Note: Given an offset and uncertainties in the photometric data set, a generalized Lomb–Scargle periodogram would be preferred (Zechmeister & Kürster 2009); however, it is not clear from the referenced papers whether this generalized method or just a basic Lomb–Scargle periodogram was used (Scargle 1982).

Notably, the estimates via a Lomb–Scargle periodogram are longer than the estimates with other methods. All results are broadly consistent with our findings from our full model results (${13.37}_{-0.17}^{+0.13}$ days; see Table 4), except the 15.2 ± 0.2 days result from the Lomb–Scargle analysis by Ciardi et al. (2018). They also have the smallest uncertainties of any rotation period estimate, so it is possible that their value is reasonable but the uncertainties are overly optimistic.

A possible explanation of this discrepancy could be differential rotation. Regardless of activity level, starspots, plage, and other activity may be more prominent at different stellar latitudes when the K2 photometry and our HARPS-N spectroscopy were conducted. This hypothesis is also mentioned by Ciardi et al. (2018) to explain a larger than expected $v\sin i$. Following Barnes et al. (2005) and Kitchatinov & Olemskoy (2012), they estimate that the equatorial rotation period of K2-136 could be faster than higher latitudes by ∼1 day. Then again, Aigrain et al. (2015) found that claims of differential rotation should be treated with caution even for long baselines of photometry. We may simply be seeing different starspots at different longitudes creating phase modulation, combined with greater or fewer numbers of starspots leading to better or worse constraints on the rotation period.

One of the planet discovery papers, Ciardi et al. (2018), reported a binary companion to K2-136. In addition to their K2 photometric analysis, they collected spectra from the SpeX spectrograph (Rayner et al. 2003, 2004) at the 3 m NASA Infrared Telescope Facility and the High Resolution Echelle Spectrometer (Vogt et al. 1994) at the Keck I telescope, as well as AO observations with the NIRC2 instrument at the Keck II telescope and the P3K AO system and the Palomar High Angular Resolution Observer camera (Hayward et al. 2001) on the 200'' Hale Telescope at Palomar Observatory. The AO observations at both facilities detected an M7/8V star separated from the primary star by ∼07, corresponding to a projected separation of ∼40 au; the spectroscopic observations did not detect this companion, and no further companions were found by any of the above observations. (Notably, this angular separation is more than enough for HST to resolve; see Sections 4.9 and 2.5.)

Gaia DR3 did not detect the binary companion, leaving the issue of boundedness unresolved. However, Ciardi et al. (2018) compared the current position of K2-136 against observations from the 1950 Palomar Observatory Sky Survey (POSS I) and noted that the star had moved 6'' in the intervening time with no evidence of background stars. These POSS I observations show that the stellar companion is likely bound.

Further, Gaia DR3 reported K2-136 to have an astrometric excess noise of 96 μas and a renormalized unit weight error (RUWE) of 1.23, a mild departure from a good single-star model. At the separation and brightness of the companion, this excess variability in the astrometry is unlikely to be due to pollution from its light contribution: at ∼07 separation a companion can be detected only with a G magnitude difference of ≲2 (Gaia Collaboration et al. 2021). There is therefore a mild indication of astrometric variability due to unmodeled orbital motion. Using the formalism of Torres (1999), at the distance of the system, given its angular separation, and for the mass range of an M7/8 star, the median astrometric acceleration is expected to be ∼25 μas yr⁻², with maximum value close to 40 μas yr⁻², indicating that in addition to simple astrometric noise, the bulk of the astrometric variability could be caused by the detection of the acceleration due to the companion.

It is worth considering whether flux from the companion could bias the measured RVs of the primary K dwarf. Both the HARPS-N and ESPRESSO bandpasses are approximately 380–690 nm and centered on the V band. According to Ciardi et al. (2018), the M dwarf companion is at least 10 mag fainter than the primary in the V band. We can use Δm = 10 as a worst-case scenario and similarly assume the companion star was well-centered on the fiber for all observations (because the companion and primary are separated by 07, the companion would not be well-centered and the actual flux contamination from the companion would be less). Cunha et al. (2013) explored the RV impact of flux contamination from a stellar companion: for a K5 dwarf and an M dwarf (M3 or later) with Δm = 10, the maximum impact on RVs is <10 cm s⁻¹ and therefore negligible for our level of RV precision.

Ciardi et al. (2018) explored whether the transit signals may originate from the M dwarf. They found that in order to match the observed transit depth of K2-136c, the M dwarf would have to be a binary system itself that exhibits significant and detectable secondary eclipses, which have not been observed. Further, the transit duration of K2-136c is inconsistent with a transit of an M dwarf. Finally, K2-136c has already been validated by Ciardi et al. (2018), and all three planets have been independently validated by M18 and Livingston et al. (2018). Therefore, it is very unlikely that the planets are false-positive signals or planetary signals from the M dwarf.

However, in the Kepler bandpass, the companion M dwarf is 6.5 mag fainter, which leads to a very small dilution effect on the planet transit depths. Following Ciardi et al. (2015), we include this effect for the sake of robustness in our final planet radius estimates (which enlarges each planet by ∼0.13%).

This binary companion may also cause a long-term RV trend, which we discuss further in Section 4.5 and test in Section 4.10.

In order to investigate periodic signals in our data (planetary or otherwise), we created generalized Lomb–Scargle periodograms (Scargle 1982; Zechmeister & Kürster 2009) of our RVs and our stellar activity indices. As a point of reference, we also included the window function of our data (built from constant, nonzero values at each of the time stamps of our observations). It is used to determine the regular patterns in the periodograms due to the sampling and gaps in the time series. These periodograms are presented in Figure 1.

Zoom In Zoom Out Reset image size

To ascertain the robustness of any apparent signals, we also estimated each periodogram's false-alarm probability (FAP), the likelihood that an apparent signal of a given strength will be detected when no underlying signal is actually present. The FAP was estimated with the bootstrap method: sampling the observations randomly with replacement while maintaining the same time stamps. We repeat this process 100,000 times, each time constructing a periodogram and determining the maximum peak. This reveals how often a given signal strength will appear due only to noise, from which the FAP is calculated.

The strongest RV signals in our combined periodogram (both HARPS-N and ESPRESSO) are at the orbital period of K2-136c, the rotation period of the star, and near 0.017 day⁻¹ (i.e., half the length of the ESPRESSO data baseline of 117.7 days), although none are significant at the 1% level. Among the stellar activity periodograms, the strongest signals are the rotation period signal in the ESPRESSO FWHM power and a long-period signal in the HARPS-N FWHM power, which can be attributed to the window function. Also, the strongest peak in the (combined data) window function periodogram above 0.01 day⁻¹ (near 0.03 day⁻¹) does not correspond to any significant signals or aliases in any of the four data types. As for low-frequency signals <0.01 day⁻¹, we discuss possible long-term trends in Sections 3.2 and 4.5. All of this indicates that K2-136c has a more detectable RV signal than the other two planets, as expected given its size.

We also examined the relationship between the RVs and our stellar activity indices. Because the RVs are measured from small shifts in spectral absorption lines and stellar activity can change the shape of absorption lines, stellar activity can significantly affect RV observations (e.g., Queloz et al. 2001; Haywood 2015; Rajpaul et al. 2015). Scatter plots between the RVs and all three activity indices are presented in Figure 2. At least in the case of S_HK and FWHM, there are notable correlations with the RVs according to the p-values for the Spearman correlation coefficient (which captures nonlinear, monotonic correlations and uses the same −1 to 1 range as the Pearson correlation coefficient). According to the Spearman coefficient, there is a correlation of 0.26 between the RV and S_HK , a negative correlation of −0.18 between the RV and BIS, and the strongest correlation of 0.39 between the RV and FWHM. For a data set of this size, these coefficients correspond to p-values of 0.004 for the RV and S_HK , 0.052 for the RV and BIS, and ≪0.001 for the RV and FWHM. Therefore, for the RV and S_HK and especially the RV and FWHM, there appears to be a statistically significant correlation. In other words, there is good reason to believe that this data set includes correlated and structured stellar activity. In fact, the Spearman correlation coefficient is likely an underestimate of the correlation between the RVs and stellar activity indices, as there can be a phase shift in the amplitude variation from one data type to another (Santos et al. 2014; Collier Cameron et al. 2019).

Zoom In Zoom Out Reset image size

We cleaned and flattened the photometric data from K2 using the exact same procedure as was used originally in M18. Their procedure follows the self-flat-fielding (SFF) method developed in Vanderburg & Johnson (2014) to perform a rough removal of instrumental variability followed by a simultaneous fit to a model consisting of Mandel & Agol (2002) transit shapes for the three planets, a basis spline in time to describe the stellar variability, and splines in Kepler's roll angle to describe the systematic photometric errors introduced by the spacecraft's unstable pointing (Vanderburg et al. 2016). After performing the fit, they removed the best-fit systematics and stellar variability components, isolating the transits for further analysis. The interested reader should refer to M18 for additional detail of the full procedure.

We modeled the flattened and cleaned M18 light curve with the BATMAN Python package (Kreidberg 2015), based on the Mandel & Agol (2002) transit model. Our model included a baseline offset parameter and white noise parameter for our K2 Campaign 13 photometry as well as two quadratic limb-darkening parameters (parameterized using Kipping 2013). Each planet was modeled with five parameters: the transit time, orbital period, planet radius relative to stellar radius, transit duration, and impact parameter. All parameters were modeled with either uniform, Gaussian, Jeffreys, or modified Jeffreys priors (Gregory 2007). Only the photometric white noise parameters used modified Jeffreys priors, with a knee located at the mean of the photometric flux uncertainty for that particular campaign or sector. All priors are listed in Table 3. The raw and flattened data can be seen in Figure 3.

Zoom In Zoom Out Reset image size

Table 3. K2-136 Transit and Planetary Parameters

Parameter	Unit	This Paper	Priors
Planet b
Period Pb	day	7.97525 ± 0.00073	Unif(7.96529,7.98529)
Time of transit t0,b	BJD-2450000	${8679.083}_{-0.074}^{+0.075}$	Unif(8678.58762,8679.58762) a
Planet–star radius ratio Rb /R*	⋯	${0.01370}_{-0.00036}^{+0.00041}$	Jeffreys(0.001,0.1)
Radius Rb	R⊕	${1.014}_{-0.049}^{+0.050}$	⋯
Transit duration T14,b	hr	${2.67}_{-0.084}^{+0.086}$	Unif(0,7.2)
Impact parameter bb	⋯	${0.22}_{-0.14}^{+0.15}$	Unif(0,1)
Semimajor axis ab	AU	0.0707 ± 0.0012	⋯
Mean density ρb	ρ⊕	<2.8b, <4.4c	⋯
Mean density ρb	g cm−3	<16b, <24c	⋯
Insolation flux Sb	S⊕	${33.5}_{-1.5}^{+1.6}$	⋯
Equilibrium temperature Teq,b (albedo = 0.3)	K	610	⋯
Equilibrium temperature Teq,b (albedo = 0.5)	K	560	⋯
Planet c
Period Pc	day	${17.30723}_{-0.00020}^{+0.00019}$	Unif(17.30514, 17.30914)
Time of transit t0,c	BJD-2450000	${8678.0792}_{-0.0096}^{+0.0088}$	Unif(8677.0747, 8679.0747) a
Planet–star radius ratio Rc /R*	⋯	${0.04064}_{-0.00071}^{+0.00068}$	Jeffreys(0.001, 0.1)
Radius Rc	R⊕	3.00 ± 0.13	...
Transit duration T14,c	hr	${3.449}_{-0.031}^{+0.039}$	Unif(0, 7.2)
Impact parameter bc	⋯	${0.31}_{-0.14}^{+0.11}$	Unif(0, 1)
Semimajor axis ac	AU	${0.1185}_{-0.0021}^{+0.0020}$	⋯
Mean density ρc	ρ⊕	${0.67}_{-0.10}^{+0.12}$	⋯
Mean density ρc	g cm−3	${3.69}_{-0.56}^{+0.67}$	⋯
Insolation flux Sc	S⊕	${11.91}_{-0.53}^{+0.57}$	⋯
Equilibrium temperature Teq,c (albedo = 0.3)	K	470	⋯
Equilibrium temperature Teq,c (albedo = 0.5)	K	440	⋯
Planet d
Period Pd	day	${25.5750}_{-0.0021}^{+0.0022}$	Unif(25.5551,25.5951)
Time of transit t0,d	BJD-2450000	${8675.936}_{-0.068}^{+0.072}$	Unif(8675.4401, 8676.4401) a
Planet–star radius ratio Rd /R*	⋯	${0.02119}_{-0.00061}^{+0.00057}$	Jeffreys(0.001, 0.1)
Radius Rd	R⊕	${1.565}_{-0.076}^{+0.077}$	⋯
Transit duration T14,d	hr	${3.04}_{-0.09}^{+0.10}$	Unif(0, 7.2)
Impact parameter bd	⋯	${0.677}_{-0.049}^{+0.042}$	Unif(0, 1)
Semimajor axis ad	AU	${0.1538}_{-0.0027}^{+0.0026}$	⋯
Mean density ρd	ρ⊕	<0.35b, <0.79c	⋯
Mean density ρd	g cm−3	<1.9b, <4.3c	⋯
Insolation flux Sd	S⊕	${7.07}_{-0.32}^{-0.34}$	⋯
Equilibrium temperature Teq,d (albedo = 0.3)	K	420	⋯
Equilibrium temperature Teq,d (albedo = 0.5)	K	380	⋯
System parameters
Kepler/K2 quadratic limb-darkening parameter q1,Kepler	⋯	${0.38}_{-0.13}^{+0.22}$	Unif(0, 1)
Kepler/K2 quadratic limb-darkening parameter q2,Kepler	⋯	${0.56}_{-0.18}^{+0.24}$	Unif(0, 1)
K2 Campaign 13 normalized baseline offset	ppm	${0.7}_{-3.0}^{+2.9}$	Unif(−1000, 1000)
K2 Campaign 13 photometric white noise amplitude	ppm	53.9 ± 1.8	ModJeffreys(1, 1000, 0)

Note.

^a t₀ is centered between the K2 and TESS photometry, so the ephemeris drift is incorporated into t₀ and P.

Download table as:  ASCIITypeset image

We also applied a Gaussian prior on stellar density by comparing the spectroscopically derived stellar density to the stellar density found via the following equation (Seager & Mallén-Ornelas 2003; Sozzetti et al. 2007):

$$\begin{eqnarray}&&{\rho }_{* }=\displaystyle \frac{3\pi }{{{GP}}^{2}}{\left(\displaystyle \frac{a}{{R}_{* }}\right)}^{3},\end{eqnarray} \tag{ 1 }$$

where the orbital period (P) and the semimajor axis (a/R_*) are derived directly from the light-curve model.

All together, our full transit model includes 15 planetary parameters (5 per planet: time of transit, orbital period, ratio of planet radius to stellar radius, transit duration, and impact parameter), 2 quadratic limb-darkening parameters, 1 photometric noise parameter, and 1 photometric baseline offset parameter for a total of 19 parameters.

The only parameters in common between our transit model and our RV model (as explained in further detail below) are the transit times and orbital periods.

No preprocessed light curves were available for the TESS Sector 43 observation of K2-136, so we extracted the photometry from the full-frame image (FFI) pixel level. Following Vanderburg et al. (2019), we constructed 20 different apertures (10 circular, 10 shaped like the TESS point-spread function) and selected the one that best minimized photometric scatter. As for TESS Sector 44 observations, we used the simple aperture photometry (SAP) light curve produced by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016). Light curves from both sectors were then flattened in the same way: a basis spline fit was performed iteratively on the photometry (with breakpoints every 0.3 day in order to adequately model stellar variability) and 3σ outliers were removed until convergence (this too, aside from the breakpoint length, follows Vanderburg et al. 2019). Finally, we conducted a simultaneous fit of the low-frequency variability and the transits in order to determine the best-fit low-frequency variability.

TESS photometry is not incorporated into our final photometric model, although we did run exploratory joint transit models on the K2 and TESS photometry simultaneously. The transit signals of the two smaller planets, K2-136b and d were too small to reliably detect in the TESS photometry: individual transits were indistinguishable in depth and quality from temporally adjacent stellar activity. However, the transits of K2-136c were easily identifiable individually in TESS photometry, so we ran a joint transit model on all K2 photometry and TESS photometry to explore the resulting improvements in the parameters of K2-136c. This joint transit model included all parameters listed in Section 4.3 as well as additional baseline offset and white noise parameters for the two TESS sectors (43 and 44) and two additional quadratic limb-darkening parameters for TESS photometry for six additional parameters. The fit resulted in consistent values for all planet and system parameters as well as a dramatically more precise ephemeris for K2-136c: ${P}_{c}={17.307081}_{-0.000013}^{+0.000014}$ days and ${t}_{0,c}={8678.07179}_{-0.00063}^{+0.00067}$ (BJD-2450000). For comparison, this period and transit time have uncertainties that are both approximately 15× tighter than those resulting from K2 transit modeling alone (see Table 3). We report these values here to minimize ephemeris drift and facilitate planning of future transit observations of K2-136c.

We modeled the RV signal of the orbiting planets and the stellar activity simultaneously. We assumed noninteracting planets with Keplerian orbits. We used RadVel (Fulton et al. 2018) to model the RV signal from each planet with five parameters: reference epoch, orbital period, RV semi-amplitude, eccentricity, and argument of periastron. The latter two parameters, the eccentricity and longitude of periastron, were parameterized as $\sqrt{e}\cos w$ and $\sqrt{e}\sin w$. As explained in Eastman et al. (2013), this reparameterization avoids a boundary condition at zero eccentricity that may lead to eccentricity estimates that are systematically biased upward. We conducted trial simulations with circular versus eccentric orbits for all three planets and found excellent agreement in all parameters (less than 1σ). We opted to keep the eccentricity and argument of the periastron as parameters in order to constrain or place upper limits on each planet's eccentricity. Additionally, we prevented system configurations that would lead to orbit crossings of any two planets, as well as overlaps of planetary Hill spheres. For each planet's reference epoch and orbital period, we applied a Gaussian prior based on the transit parameters determined from M18. We analyzed the K2 transit photometry with and without TESS photometry and verified these M18 values (see Table 3).

We also applied additional prior limits on the eccentricities of K2-136b and K2-136d. Based on preliminary modeling of all three planets with uninformed prior eccentricity constraints (and everything else identical to our final model), we found that we could determine the eccentricity of K2-136c but not its siblings. Thus, we decided to set eccentricity constraints by using the Stability of Planetary Orbital Configurations Klassifier (SPOCK; Tamayo et al. 2020), an N-body simulator that employs machine learning to improve performance. We input into SPOCK our stellar mass posterior and our orbital period posteriors for all three planets (from our preliminary simulations). We also input the mass, eccentricity, and argument of the periastron posteriors for K2-136c (again, from our early simulations). Lastly, we input uniform distributions for the mass, eccentricity, and argument of the periastron for K2-136b and K2-136d. The planet mass ranged from 0 up to that of approximately 100% iron planet composition (3 M_⊕ for K2-136b and 30 M_⊕ for K2-136d; Fortney et al. 2007). The eccentricity ranged from 0 to 0.6 and argument of the periastron ranged from 0° to 360°. We took the subset of our sample that had a >90% chance of surviving for >10⁹ orbits (of the innermost planet) and then determined the 3σ upper limit on the eccentricities of K2-136b and K2-136d for that subset. We then used those values, ${e}_{{\rm{b}},\max }=0.35$ and ${e}_{{\rm{d}},\max }=0.37$, as eccentricity upper limit priors in all subsequent simulations.

Given the M dwarf companion to our host star (Ciardi et al. 2018), we wanted to include the potential for an RV trend caused by this companion. For further discussion of binarity and linear trends, see Section 3.2. Our planetary RV signals and the RV trend therefore take the following form:

$$\begin{eqnarray}&&{RV}=\displaystyle \sum _{i\in \{b,c,d\}}{K}_{i}\left(\cos \left({\omega }_{i}+{f}_{i}\right)+{e}_{i}\cos {\omega }_{i}\right)+{mt},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&f=2\arctan \left(\sqrt{\displaystyle \frac{1+e}{1-e}}\tan \displaystyle \frac{E}{2}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&M=E-e\sin E,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&M={nt}=\displaystyle \frac{2\pi (t-\tau )}{P},\end{eqnarray} \tag{ 5 }$$

where K is the induced RV semi-amplitude, ω is the argument of the periastron, f is the true anomaly, e is the eccentricity, m is the slope of the RV trend, t is the observation time, E is the eccentric anomaly, M is the mean anomaly, n is the mean motion, τ is the time of periastron passage (as calculated from transit time, orbital period, eccentricity, and argument of the periastron), and P is the orbital period.

The stellar activity was handled via a GP on the RVs. We first used a simultaneous model of stellar activity on four different data types: RV, FWHM (a measure of the width of absorption lines), BIS (a measure of line asymmetry), and S_HK (an estimate of chromospheric magnetic activity via emission in the cores of the Ca ii H & K lines). However, we found that this approach forced the model to include unreasonable amounts of white noise into each data type via a white noise jitter parameter included in our model. Further, it did not lead to a notable improvement in our final parameter constraints. We conducted numerous tests exploring the excess white noise preferred by the model, including modeling different instruments, different numbers of planets, altering the data reduction pipeline (e.g., trying the ZLSD pipeline; Lienhard et al. 2022), and creating synthetic data sets to model against. The most reasonable hypothesis we could find is that for K2-136, our data (especially for ESPRESSO) are of a sufficiently high quality, with small enough uncertainties, that the model we used from Rajpaul et al. (2015) to relate the stellar activity indices to each other and to the RVs was not complex enough to account for the correlated structure of the stellar activity of K2-136.

Our RV data and the final RV model fit can be seen in Figure 4. Despite only modeling RVs without any stellar activity indices, a GP is still robust and allows us to separate planet-induced RVs from stellar activity to the extent that disentanglement is possible. We followed the method laid out in Rajpaul et al. (2015), but we constrained the model only to the portions relevant to RVs rather than additional stellar activity indices. They dictate that the RVs are related to the stellar activity as follows:

$$\begin{eqnarray}&&{\rm{\Delta }}{RV}={V}_{c}G(t)+{V}_{r}\dot{G}(t).\end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

In this equation, G(t) corresponds to the underlying stellar activity GP, while V_c and V_r correspond to the RV amplitudes of the convective blueshift and rotation modulation effect, respectively. It is important to include both the rotation modulation and convective blueshift for data types impacted by both, as one phenomenon may play a larger role than the other depending on the data type and star. In fact, our final results (see Table 4) show that for K2-136, rotation modulation has an outsized effect on RVs compared to convective blueshift; however, we retain both terms in our model as V_c is still inconsistent with zero. We follow Rajpaul et al. (2015) further by using a quasiperiodic kernel to establish the covariance matrix of our GP. A quasiperiodic kernel is a good choice to capture the stellar variability of a star because quasiperiodicity describes well the variability exhibited by a rotating star with starspots that come and go. Following Giles et al. (2017), we find a predicted starspot lifetime for K2-136 of ${38}_{-13}^{+20}$ days (versus a rotation period of 13.37 days); this estimate is consistent with our GP evolution timescale result (${48}_{-11}^{+19}$ days, see Table 4). In other words, stellar activity during the time frame of a single rotation period is likely to look similar (though not identical) to stellar activity during the previous and subsequent stellar rotation periods, as the starspots for K2-136 are likely to be longer lived than the stellar rotation period. Thus, the quasiperiodic kernel is defined as follows:

$$\begin{eqnarray}&&\begin{array}{l}{K}_{\mathrm{QP}}({t}_{i},{t}_{j})\\ =\,{h}^{2}\exp \left(-\displaystyle \frac{{\sin }^{2}\left(\pi ({t}_{i}-{t}_{j})/{P}_{* }\right)}{2{\lambda }_{p}^{2}}-\displaystyle \frac{{\left({t}_{i}-{t}_{j}\right)}^{2}}{2{\lambda }_{e}^{2}}\right),\end{array}\end{eqnarray} \tag{ 7 }$$

where h is the GP amplitude (which is folded into the amplitude parameters described above), P_* is the stellar stellar rotation period, λ_e is a decay timescale proportional to the starspot lifetime, and λ_p is a smoothness parameter that captures the level of variability within a single rotation period. t_i and t_j are any two times between which the covariance is calculated; for a given time series of N observations, all N² combinations of time pairs create the NxN covariance matrix. This covariance matrix (plus a mean model) is a normal multivariate distribution; G(t) can be explored by sampling from this distribution. We refer the reader to Rajpaul et al. (2015) for a more detailed description of this method, as well as Mayo et al. (2019) for an application of this method to the sub-Neptune Kepler-538b.

Table 4. K2-136 RV Model Parameters

Parameter	Unit	This Paper	Priors
Planet b
Period Pb	day	7.97520 ± 0.00079	Normal(7.97529, 0.00080) a
Time of transit t0,b	BJD-2450000	${7817.7563}_{-0.0048}^{+0.0046}$	Normal(7817.7563, 0.0048) a
Semi-amplitude Kb	m s−1	<1.2 b , <1.7 c	Unif(0.001, 20)
Eccentricity eb	⋯	${0.14}_{-0.11}^{+0.12}$ (<0.21 b , <0.32 c )	d d , e e
Argument of periastron ωb	degrees	${189}_{-132}^{+82}$	d, e
Mass Mb	M⊕	<2.9 b , <4.3 c	⋯
Planet c
Period Pc	day	17.30713 ± 0.00027	Normal(17.30714, 0.00027) a
Time of transit t0,c	BJD-2450000	${7812.71770}_{-0.00085}^{+0.00086}$	Normal(7812.71770, 0.00089) a
Semi-amplitude Kc	m s−1	${5.49}_{-0.52}^{+0.54}$	Unif(0.001, 20)
Eccentricity ec	⋯	${0.047}_{-0.034}^{+0.062}$ (<0.074 b , <0.16 c )	d
Argument of periastron ωc	degrees	124 ± 99	d
Mass Mc	M⊕	${18.1}_{-1.8}^{+1.9}$	⋯
Planet d
Period Pd	day	${25.5750}_{-0.0023}^{+0.0024}$	Normal(25.5751, 0.0024) a
Time of transit t0,d	BJD-2450000	7780.8117 ± 0.0065	Normal(7780.8116, 0.0065) a
Semi-amplitude Kd	m s−1	<0.36 b , <0.78 c	Unif(0.001, 20)
Eccentricity ed	⋯	${0.071}_{-0.049}^{+0.063}$ (<0.10 b , <0.16 c )	d, e
Argument of periastron ωd	degrees	${280}_{-110}^{+130}$	d, e
Mass Md	M⊕	<1.3 b , <3.0 c	⋯
System and GP parameters
RV slope	m s−1 yr−1	0.1 ± 2.5	Unif(−365, 365)
HARPS-N RV white noise amplitude	m s−1	0.83 ± 0.52	Unif(0, 20)
ESPRESSO RV white noise amplitude	m s−1	${1.57}_{-0.62}^{+0.73}$	Unif(0, 20)
HARPS-N RV offset amplitude	m s−1	${39503.7}_{-1.9}^{+2.1}$	Unif(39450, 39550)
ESPRESSO RV offset amplitude	m s−1	${39494.7}_{-3.6}^{+3.5}$	Unif(39450, 39550)
GP RV convective blueshift amplitude Vc	m s−1	${3.5}_{-1.2}^{+2.2}$	Unif(0, 100)
GP RV rotation modulation amplitude Vr	m s−1	${22.0}_{-8.0}^{+12.9}$	Unif(0, 100)
GP stellar rotation period P*	day	${13.37}_{-0.17}^{+0.13}$	Unif(1, 20)
GP inverse harmonic complexity λp	⋯	${0.75}_{-0.16}^{+0.23}$	Unif(0.1, 3)
GP evolution timescale λe	day	${48}_{-11}^{+19}$	Unif(1, 200)

Notes.

^a Mann et al. (2018). ^b 68% confidence limit. ^c 95% confidence limit. ^d Unif(−1, 1) on $\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$. See Eastman et al. (2013). ^e K2-136b and K2-136d had additional eccentricity prior upper limits of 0.35 and 0.37, respectively; see Section 4.5.

Download table as:  ASCIITypeset image

Finally, we include a jitter parameter (added in quadrature to RV uncertainties) and a baseline offset parameter for each instrument (HARPS-N and ESPRESSO). All together, our full model includes 21 planetary parameters (5 per planet: time of transit, orbital period, RV semi-amplitude, $\sqrt{e}\sin \omega$, and $\sqrt{e}\cos \omega$), 5 GP parameters (2 GP amplitudes corresponding to V_c and V_r from Equation (6) as well as P_*, λ_p , and λ_e from Equation (7)), 1 linear trend parameter, 2 RV noise parameters (1 per instrument), and 2 RV baseline offset parameters (1 per instrument), for a total of 25 parameters.

We used a uniform prior for the RV semi-amplitude for each planet. However, we conducted a trial simulation to compare a uniform versus log uniform prior for the RV semi-amplitude and found no discernible difference in our results (all parameters agreed to within 1σ).

We conducted parameter estimation of our model with the observed data using the Bayesian inference tool MultiNest (Feroz et al. 2009, 2019). We set MultiNest to the constant efficiency mode, importance nested sampling mode, and multimodal mode. We used a sampling efficiency of 0.01, 1000 live points, and an evidence tolerance of 0.1. Constant efficiency is typically off, but we turned it on as it allows for better exploration of parameter space in higher-dimensional models such as our own. Further, sampling efficiency is usually set to 0.8, and the number of live points is usually set to 400. Decreasing sampling efficiency and increasing the number of live points leads to more complete coverage of parameter space, at the cost of a typically longer simulation convergence time. Finally, the evidence tolerance is usually set to 0.8; reducing the evidence tolerance causes the simulation to run longer but increases confidence that the simulation has fully converged. In other words, the standard MultiNest settings would likely lead to reliable results, but our choice of settings increases the trustworthiness of our parameter estimation and model evidence results.

One of the strengths of MultiNest is that it automatically calculates the Bayesian evidence of the selected model, making model comparison very easy. We compared the model evidences of eight different models (RV only, no photometry), based on all possible combinations of planets b, c, and d. The results of our model comparisons are listed in Table 5. We find that the most preferred model is the one that contains only planet c, and not planets b or d. In fact, using the Bayes factor interpretation of Kass & Raftery (1995), we find that almost every other combinations of planets can be decisively ruled out (i.e., the Bayes factor of the planet c model to any other model is >100). The only exceptions are the model with planets b and c and the model with planets c and d, which are only strongly disfavored. This tells us that neither K2-136b nor K2-136d are unambiguously detected, so including either in the model quickly worsens the model evidence. However, it is notable that the model with planets b and c is better than the model with planets c and d, and is nearly in the more likely "Disfavored" category rather than "Strongly disfavored" one. This makes sense, as K2-136b (unlike K2-136d) has a nonzero peak in its semi-amplitude posterior distribution (see Figure 5). In fact, although only upper limits are reported for the mass of K2-136b in Table 4, the median mass is actually nonzero at the 2.0σ level (M_b = 2.4 ± 1.2 M_⊕).

Zoom In Zoom Out Reset image size

Table 5. Model Evidence Comparisons for K2-136 Planet Configurations

Planets in Model	${\rm{\Delta }}{\mathrm{log}}_{10}$(Evidence)	Interpretation
c	0	⋯
b,c	−1.10	Strongly disfavored
c,d	−2.06	Strongly disfavored
b,c,d	−3.23	Excluded
⋯	−5.54	Excluded
d	−7.00	Excluded
b	−7.78	Excluded
b,d	−9.17	Excluded

Download table as:  ASCIITypeset image

Our results tell us that we can be confident that K2-136c has been detected in our observations. In contrast, the RV signals of K2-136b and K2-136d fall below the threshold of detection, at least given the quantity and quality of our specific data set (Figure 6 still includes planetary mass and density upper limits for K2-136b and K2-136d). Continued radial velocity monitoring, particularly with high-precision instruments and facilities may be able to measure their masses, especially K2-136b (which appears to already be near the threshold of detection with our current data set).

Although the model with only K2-136c is the most favored, we still use a model with all three planets as our canonical model for parameter estimation for three reasons. First, we already know from transit photometry that planets b and d exist. Accordingly, the goal of the model comparison exercise described above is not to question the existence of these planets but to examine whether the RV signals from each planet can be detected in our data set. By adopting the three-planet model, we incorporate the uncertainties introduced by the unknown masses and eccentricities of planets b and d. Second, using the three-planet model allows us to determine upper limits for the masses of planets b and d, which is useful for constraining planet compositions and providing guidance for any future attempts to constrain the mass of either planet. Third, both models agree very closely: all parameters are consistent at 1σ or less, and all uncertainties on parameters are similar in scale. For the RV semi-amplitude of K2-136c, our key parameter of interest, our canonical model returned ${K}_{c}={5.49}_{-0.52}^{+0.54}$ m s⁻¹ while the one-planet model returned ${K}_{c}={5.17}_{-0.51}^{+0.56}$ m s⁻¹.

In order to rigorously assess the accuracy of our results, we conducted tests to analyze different components of our model. In particular, we removed the GP portion of the model, and we also injected and recovered synthetic planets into the system to compare input and output RV semi-amplitudes.

For our test models we chose not to include planets b and d, as well as photometry, and we then compared against the model with only K2-136c (hereafter referred to as the "one-planet reference model") rather than the three-planet model; this is despite already selecting the three-planet model as the canonical model to report our results (see Section 4.7). This was done for a few reasons. First, the one-planet model is the preferred model according to the Bayesian evidences, so it is a reasonable point of comparison. Second, as stated earlier, the one-planet reference model and the three-planet model agree very closely. Therefore, any test model parameters found to be highly consistent with the one-planet reference model results will also be highly consistent with the three-planet model results. Third, as a practical matter, including only K2-136c in our test models significantly reduced computational complexity, allowing us to test a wider variety of models.

4.8.1. No GP

4.8.2. Synthetic Planet Injections

4.9.1. HST NUV Observations

To compare the UV quiescent activity of K2-136 with other K stars Hyades members, we measured the surface flux of the Mg ii h (2796.35 Å) and k (2808.53 Å) lines. The Mg ii lines are the strongest emission lines in the NUV and correlate strongly with the chromospheric activity of the star. For accurate emission measurements, we subtracted the NUV continuum by fitting the data outside of the Mg ii integration region using the astropy module specutils. We then integrated over 2792.0–2807.0 Å to measure the Mg ii emission flux. To convert from observed flux to surface flux, we estimated the radii of the stars using the relationships between age, effective temperature, and radius given by Baraffe et al. (2015), except for K2-136, where we use the radius reported in this work.

Zoom In Zoom Out Reset image size

Figure 7 shows the Mg ii surface flux as a function of the rotation period for K2-136 and 13 other observed K star members of the Hyades from GO-15091. K2-136 has a surface Mg ii flux of (8.32 ± 0.17) × 10⁵ erg s⁻¹ cm⁻², whereas the median of the sample is (9.66 ± 0.24) × 10⁵ erg s⁻¹ cm⁻².

Zoom In Zoom Out Reset image size

Additionally, Richey-Yowell et al. (2019) measured the NUV flux densities of 97 K stars in the Hyades using archival data from the Galaxy Evolution Explorer (GALEX; Martin et al. 2005). The median NUV flux density of their sample of Hyades stars at a normalized distance of 10 pc was 1.89 × 10³ μ Jy. We measure a GALEX NUV magnitude for K2-136 of 19.51 mag, which corresponds to a flux density at 10 pc of 1.52 × 10³ μ Jy, well within the interquartiles of the total sample from Richey-Yowell et al. (2019).

4.9.2. X-Ray Observations

We detected a total of 800 European Photon Imaging Camera (EPIC) counts for K2-136 in the XMM observation and can therefore extract an X-ray spectrum for the star. We used a one-temperature Astrophysical Plasma Emission Code (APEC) model to fit the spectrum, which is appropriate for representing the hot plasma in stellar coronae. We combined this model with the interstellar medium absorption model tbabs using photoelectric cross section from Balucinska-Church & McCammon (1992) to account for the neutral hydrogen column density N_H. We set N_H to 5.5 × 10¹⁸ cm⁻², derived using E(B − V) = 0.001 for Hyades (Taylor 2006), R_V = 3.1, and the relation N_H[cm⁻²/A_v ] = 1.79 × 10²¹ (Predehl & Schmitt 1995); allowing N_H to float did not improve the fit. The spectra for each of the three XMM cameras are shown in Figure 8 and the best-fit parameters, which are obtained from a simultaneous fit to the three, are provided in Table 6.

Zoom In Zoom Out Reset image size

Table 6. X-Ray Spectral Fit Parameters of K2-136

Parameter a	Value	Unit
Degrees of Freedom	33
Reduced χ2	0.91
Plasma Temperature	0.65 ± 0.07	keV
Plasma Metal Abundance	0.06 ± 0.02	Solar
pn Energy Flux	3.12 ± 0.68	10−14 erg s−1 cm−2
MOS1 Energy Flux	3.01 ± 0.93	10−14 erg s−1 cm−2
MOS2 Energy Flux	3.02 ± 0.88	10−14 erg s−1 cm−2
EPIC b Energy Flux	3.06 ± 0.46	10−14 erg s−1 cm−2

Notes.

^a All flux values are in the 0.1−2.4 keV energy range. ^b The error-weighted average of the pn and MOS cameras.

Download table as:  ASCIITypeset image

Using the total EPIC energy flux from the spectral best fit, we obtained an X-ray luminosity L_X = (1.26 ± 0.19) × 10²⁸ erg s⁻¹ and L_X/L_* = (1.97 ± 0.30) × 10⁻⁵ (0.1–2.4 keV energy range). These values are within 1σ of those found by Fernandez Fernandez & Wheatley (2022) using the same XMM observation. ³⁰ For the sample of 89 K dwarfs with X-ray detections in the Hyades, the median values for L_X and L_X/L_* are ${4.5}_{-3.3}^{+7.9}\times {10}^{28}$ erg s⁻¹ and ${4.0}_{-1.6}^{+20.8}\times {10}^{-5}$, respectively (Núñez et al. 2022). K2-136, therefore, appears somewhat less luminous in X-rays than most of its coeval K dwarf brethren. A narrower comparison, against late-K (K5 and later) dwarf Hyads, shows that the L_X and L_X/L_* values for K2-136 are within one standard deviation of the median for that cohort.

In addition to X-ray luminosity, we also estimated extreme ultraviolet (EUV) luminosity using stellar age and Equation (4) from Sanz-Forcada et al. (2011) and found L_EUV = (${22.6}_{-5.7}^{+7.8}$) × 10²⁸ erg s⁻¹. Combining L_X and L_EUV and using the semimajor axis of K2-136c, we are able to estimate the X-ray and UV flux incident on K2-136c to be F_XUV = (${6.0}_{-1.4}^{+2.0}$) × 10³ erg s⁻¹ cm⁻². Then, using this incident flux value (along with M_c and R_c) we estimate an atmospheric mass-loss rate with Equation (1) from Foster et al. (2022). This yields a current atmospheric mass-loss rate for K2-136c of $\dot{{M}_{{\rm{c}}}}$ = (${3.4}_{-0.9}^{+1.3}$) × 10⁹ g s⁻¹ = (${17.8}_{-5.0}^{+6.7}$) × 10⁻³ M_⊕ Gyr⁻¹. This rate is based on current values, so the mass-loss rate in the past or future may differ. At this current rate, with a H₂–He envelope mass fraction of ∼5%, it would take ${51}_{-16}^{+23}$ Gyr to fully evaporate the atmosphere. Even the 95% lower limit evaporation time is still 28 Gyr, longer than the age of the universe. In other words, in ∼4 Gyr, when the K2-136 system is as old as the solar system currently is, we expect K2-136c will have likely only lost 5%–10% of its current atmosphere.

We also calculated the Rossby number R_o of K2-136, which is defined as the star's rotation period P_* divided by the convective turnover time τ. We used the (V − K_s )-log τ empirical relation in Wright et al. 2018 (their Equation (5)) to obtain τ = 22.3 days for K2-136. Using our measured P_* value (see Section 3.1) gives R_o = 0.6. This R_o puts K2-136 well within the X-ray unsaturated regime, in which the level of magnetic activity decays follows a power slope as a function of R_o (see Figure 3 in Wright et al. 2018). For the sample of 51 K dwarf rotators with X-ray detection in the Hyades, the median R_o = ${0.46}_{-0.08}^{+0.06}$ (Núñez et al 2023, in preparation), which suggests that the lower levels of X-ray emission from K2-136 relative to its fellow Hyades K dwarfs can be explained by its slower rotation rate.

In conclusion, K2-136 does not appear unusually active in either the NUV or the X-ray relative to its fellow Hyades K dwarfs.

As discussed in Section 3.2, prior observations of this system revealed a likely bound M7/8V companion with a projected separation of approximately 40 au. The presence of a nearby stellar companion can easily cause a trend in RV observations. We wanted to estimate the range of possible trend amplitudes for our system using some reasonable assumptions. Using the distance of ${58.752}_{-0.072}^{+0.061}$ pc from Bailer-Jones et al. (2021) and the projected angular separation of 0730 ± 0030 from adding in quadrature the R.A. and decl. separation components in Ciardi et al. (2018), we found a projected separation of 42.9 ± 1.7 au. We considered the possibility of additional radial separation by folding in a uniform distribution on radial separation between 0 and twice the median projected separation to estimate an overall separation (this broad range was chosen to include radial separations of approximately the same scale as the projected separation). As an approximation, we treat this overall separation as the semimajor axis.

Next, the stellar companion was reported in Ciardi et al. (2018) to have a spectral type consistent with M7/8 (we were unable to find uncertainties associated with this result, but nearby spectral types were never mentioned). Taking a conservative approach, we assumed a stellar companion mass between 1 M_Jup (for a low-mass brown dwarf) and 0.2 M_⊙ (for a mid to late M dwarf).

Then, we combined our separation distribution and stellar mass distributions (primary and companion) via Kepler's third law to get a broad orbital period estimate of ${520}_{-180}^{+320}$ yr. Including a wide range of eccentricities (Unif(0,0.9)), we estimated an RV semi-amplitude of ${490}_{-320}^{+340}$ m s⁻¹. On the timescale of our observations, this centuries-long sinusoidal signal would manifest as a linear trend, with a maximum (absolute) slope of ${5.3}_{-3.6}^{+7.5}$ m s⁻¹ yr⁻¹. This is a very rough estimate with many assumptions, but it serves to demonstrate that a drift of only a few m s⁻¹ each year or less is very reasonable. Indeed, from our model of the RV data we found an RV trend of 0.1 ± 2.5 m s⁻¹ yr⁻¹, highly consistent with both a zero trend as well as our estimate calculated here.

To explore the suitability of the K2-136 system planets for atmospheric characterization, we calculated the transmission spectroscopy metric (TSM) defined in Kempton et al. (2018). K2-136c has a TSM of ${32.7}_{-4.1}^{+4.8}$, which is well below the recommended TSM of 90 for R_p > 1.5 R_⊕. Because K2-136b and K2-136d have unconstrained masses, we followed Zeng et al. (2016) and assumed an Earth-like CMF of 0.325 in order to predict planet masses of ${0.85}_{-0.21}^{+0.27}$ M_⊕ and ${3.35}_{-0.84}^{+1.08}$ M_⊕, respectively. For K2-136b, this yields a TSM of ${4.20}_{-0.37}^{+0.42}$, well below the recommended TSM of 10 for R_p < 1.5 R_⊕. As for K2-136d, its radius of R_d = 1.565 ± 0.077 R_⊕ is very near 1.5 R_⊕, where the TSM metric includes a scale factor that jumps dramatically, thus creating a TSM bimodal distribution. Thus, for R_d < 1.5 R_⊕ we find a TSM of ${2.36}_{-0.13}^{+0.15}$ and for R_d > 1.5 R_⊕ we find a TSM of ${13.6}_{-1.1}^{+1.0}$. In their respective radius ranges, these values are both well below the recommended TSM, so K2-136d is also probably not a good target for atmospheric characterization.

We also calculated the emission spectroscopy metric (ESM) for K2-136b and K2-136d as defined in Kempton et al. (2018). The metric applies only to "terrestrial" planets with R_p < 1.5 R_⊕, excluding K2-136c. We find K2-136b and K2-136d have ESM metrics of ${0.520}_{-0.047}^{+0.049}$ and ${0.342}_{-0.041}^{+0.043}$, respectively, both below the recommended ESM of 7.5 or higher. Therefore, these planets do not appear to be particularly attractive targets for emission spectroscopy or phase curve detection.

The TSM analysis of K2-136c is not very favorable for atmospheric characterization, but we decided to conduct a more thorough transmission spectroscopy analysis. We used the JET tool (Fortenbach & Dressing 2020) to model atmospheric spectra and to simulate the performance of the JWST instruments for certain atmospheric scenarios. We opted for the broad wavelength coverage of combining Near Infrared Imager and Slitless Spectrograph (NIRISS) Single Object Slitless Spectroscopy (SOSS) Order 1 (0.81–2.81 μm) with NIRSpec G395M (2.87–5.18 μm), as recommended by Batalha & Line (2017) to maximize the spectral information content. A single instrument, the Near Infrared Spectrograph (NIRSpec) Prism, can also cover this wavelength range, but brightness limits preclude its use here. We assumed pessimistic prelaunch instrumental noise values (Rigby et al. 2023), so a future JWST program for atmospheric characterization should outperform our conservative expectations.

The JET tool found that for an optimistic, cloudless, low-metallicity atmosphere (5× solar) we can meet a ΔBIC detection threshold of 10 (corresponding to a ∼3.6σ detection of the atmosphere compared to a flat line) with five free retrieval parameters (i.e., recon level) with only one transit for NIRISS SOSS and two transits for NIRSpec G395M. For a less optimistic, higher-metallicity atmosphere (100× solar) with clouds at 100 mbar, we can meet the same detection threshold with two transits for NIRISS SOSS and five transits for NIRSpec G395M.

With 10 free retrieval parameters (a more typical number), and with the optimistic atmosphere, we can meet a ΔBIC detection threshold of 10 with only one transit for NIRISS SOSS and three transits for NIRSpec G395M. For the less optimistic atmosphere, we can meet the detection threshold with three transits for NIRISS SOSS, but we show no detection for NIRSpec G395M for up to 50 transits considered. This analysis makes the conservative assumption that the instrument noise floor is not reduced by co-adding transits.

The resulting spectra for the low-metallicity, cloudless, case are shown in Figure 10. It seems that atmospheric characterization of K2-136c may be within reach (assuming a relatively low mean molecular weight/low-metallicity atmosphere, and low cloud level), but could require a more significant investment of JWST resources if the actual atmospheric properties are less favorable.

Zoom In Zoom Out Reset image size

It should be noted that given the on-sky position of K2-136, the ability to observe the system with JWST will be limited due to aperture position angle constraints. In addition, the very close (∼07) stellar companion may cause contamination of spectra from both instruments. This is a common issue for NIRISS as it is slitless, but NIRSpec can usually isolate the primary target with its 16 square aperture. For K2-136 the companion star is well inside this aperture boundary and will likely create some contamination. The companion is significantly fainter than the host star (J magnitude of 14.1 versus 9.1), which should mitigate the impact to a degree. It should also be possible to reduce the companion M-star spectral contamination effect post-processing.

The most enticing feature of K2-136 is the young age of the system; it could be argued that despite the potential difficulty in observing the planets' atmospheres, the rewards outweigh the risks for the chance to better understand the atmospheres of very young, relatively small planets. Observations of this system could help us construct a picture of the environment and evolution of young, low-mass planets.